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Abstract. We present similarity solutions for adiabatic bubbles that are blown by winds having time independent 
mechanical luminosities and that are each mass-loaded by the hydrodynamic ablation of distributed clumps. The 
mass loading is 'switched-on' at a specified radius (with free-expansion of the wind interior to this point) and 
injects mass at a rate per unit volume proportional to M r x where 5 = 4/3 (1) if the flow is subsonic (su perso nic 



with respect to the clumps. In the limit of negligible mass loading a similarity solution found by Dyson ( 1975 ) for 
expansion into a smooth ambient medium is recovered. The presence of mass loading heats the flow, which leads 
to a reduction in the Mach number of the supersonic freely-expanding flow, and weaker jump conditions across 
the inner shock. In solutions with large mass loading, it is possible for the wind to connect directly to the contact 
discontinuity without first passing through an inner shock, in agreement with previous hydrodynamic simulations. 
In such circumstances, the flow may or may not remain continuously supersonic with respect to the clumps. For 
a solution that gives the mass of swept-up ambient gas to be less than the sum of the masses of the wind and 
ablated material, A < —2, meaning that the exponent of the density profile of the interclump medium must be 
at most slightly positive, with negative values preferred. Maximum possible values for the ratio of ablated mass 
to wind mass occur when mass loading starts very close to the bubble center and when the flow is supersonic 
with respect to the clumps over the entire bubble radius. Whilst mass loading always reduces the temperature of 
the shocked wind, it also tends to reduce the emissivity in the interior of the bubble relative to its limb, whilst 
simultaneously increasing the central temperature relative to the limb temperature. The maximum temperature 
in the bubble often occurs near the onset of mass loading, and in some cases can be many times greater than 
the post-inner-shock temperature. Our solutions are potentially relevant to a wide range of astrophysical objects, 
including stellar wind-blown bubbles, galactic winds, starburst galaxy superwinds, and the i mpact of an AGN 
wind on its surrounding environment. This work complements the earlier work of Pittard et al. (2001) in which it 
was assumed that clumps were evaporated through conductive energy transport. 
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1. Introduction 

Numerous observations of wind-blown bubbles (WBBs) 
have led to the conclusion that their structures and evo- 
lution are significantly affected by mass entrainment from 
embedded clumps. As a case in point, bubbles blown by 
the winds of highly evolved stars can be mass-loaded from 
clumps of stellar material ejected during prior stages of 
mass loss (S mith et al. 1984 ; Har tquist et al. 1986 ; Meaburn 
et al . 1991 , Dyson & Hartquist 1992| ; Hartquist & Dyson 
1993| ). The properties of winds mass- loaded by material 
from clumps (or stellar sources) has also received theo- 
retical attention in the context of ultracompact H II re- 



gions (Dyson et al. 1995), globular cluster winds (Durisen 
& Burns 1981), galactic winds and starburst galaxy su- 
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perwinds (Strickland & Stevens 2000), and accretion flow 
structures at the centers of active galactic nuclei (David 
& Durisen |l989| ; Toniazzo et al. |2001fl . 

Investigations into the self-similar nature of super- 
nova remnant (McKee & Ostriker 1977; Chieze & Lazareff 
1981; Dyson & Hartquist |1987) and stellar wind-blown 



bubble (Pittard et al. 2001) evolution in tenuous media 
with embedded clumps have been performed. Detailed 
one-dimensional, time-dependent hydrodynamic models of 
specific WBBs associated with evolved stars have also 
been constructed (Arthur et al.|1993|, [1994 |1996|). In 



at least some isothermal mas s-loaded w inds, Arthur et 
al. (|1994|) and Williams et al. ( |l995i |l999|) showed that a 
global shock does not occur in the region where the Mach 
number is large. Instead, a global shock occurs only after 
the wind has travelled far enough for mass loading to de- 
celerate it sufficiently that its Mach number is less than 
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a few. In some of these cases, a global shock within the 
mass loading region does not exist at all, and the wind 
continues until it encounters a termination shock caused 
by the interaction of the wind with an external medium. 
Observations supporting mass loading in WBB include 
the detection of blue-shifted absorption features of species 
with a range of ionization potentials in the spectrum of the 
central star of the Wolf-Rayet nebula RCW 58. Smith et 
al. (1984) argued that the observed velocity spread for the 
detected features is much larger than would be expected 
for T < 10 5 K gas consisting only of stellar wind material 
decelerated through a terminal shock, and suggested that 
the velocity spread originates due to mass loading of the 
shocked wind through the entrainment of material from 
clumps. Crucially, one can obtain an estimate for the ratio 
of wind mass to entrained clump mass: the observed ve- 
locity spread implies a value of 40 or 50 to one. A detailed 
time-dependent, hydrodynamic model including nonequi- 
librium ionization structure is consistent with these con- 



clusions (Arthur et al. 1996 ). Spectroscopic data also pro- 
vide evidence for transonic flows in the halo of core-halo 
planetary nebula (Meaburn et al. |l99l| ), again consistent 
with high mass loading rates. Finally, it is also suspected 
that the gradual deceleration of a wind by mass loading 
and the associated weakening of an inner shock may in fact 
contribute to the radio quietness of some WBBs, although 
this is a conjecture which currently cannot be proven (see 
Williams et al. |1995| |1999 [). 

In Pittard et al. (2001, hereafter PDH), similarity so- 
lutions were derived for the structures and evolution of 
mass-loaded WBBs, under the assumption that conduc- 
tive evaporation from embedded clumps was the dom- 
inant mass loading process, and that both evaporation 
from the cold swept-up shell and radiative losses were 
negligible. To obtain a similarity solution with these con- 
ditions, specific radial power-laws on the clump and in- 
tcrclump density distribution, and temporal power-laws 
on the wind mass-loss rate and terminal velocity, were 
required. Approximate similarity solutions for evapora- 
tively mass-loaded WBBs with the assumptions of con- 
stant mass-loss rate and wind velocity, and an isobaric 
shocked wind region have previously been obtained by 
Weaver et al. ( 1977 ; evaporation from the cool swept-up 
shell) and Hanami & Sakashita (1987; mass- loading from 
clumps). A central assumption in both of these papers 
was that the shocked wind was approximately isobaric. 
However, PDH demonstrated that this was not necessarily 
a good assumption (e.g. see their Fig. 4). Indeed, impos- 
ing this condition is likely to set a limit to the amount of 
mass loading that can occur. 

A central conclusion of PDH was that there exist max- 
imum possible values for the ratio of evaporated mass to 
stellar wind mass, as a consequence of the evaporation 
rates dependence on temperature and the lowering of tem- 
perature by mass loading. In particular it was difficult to 
find ratios approaching what was observationally required. 
The work in this paper complements PDH by considering 
the case in which hydrodynamic ablation is the dominant 



mass addition process. As conductively driven evaporation 
has a very temperature sensitive rate, ablation is likely to 
regulate clump dispersal into lower temperature media. 

Our solutions are again potentially relevant to such 
diverse objects as WBBs created by a faster wind inter- 
acting with a clumpy AGB superwind, by the wind of a 
young star interacting with surrounding molecular mate- 
rial, and the wind of an active galactic nucleus impacting 
its environment. 

2. Similarity Solutions 

The basic physics behind the hydrodynamic ablation of 
material from dense clumps i nto th e surrounding flow was 
presented by Hartquist et al. ( |l986|) . It was proposed that 
if the flow around a clump is subsonic, hydrodynamic mix- 
ing occurs as a result of the well-known Bernoulli effect. 
For flow with a Mach number M with respect to the clump 
of less than unity, this leads to a volume mass injection 
rate, p, proportional to M 4 ' 3 . For supersonic flow, mixing 
occurs largely as a result of a low pressure region over the 
reverse face of the clump, 'shadowed' from the flow. In this 
case the mass injection rate is Mach number independent. 

In this paper we present similarity solutions for WBBs 
mass-loaded by hydrodynamic ablation. Our solution 
method is very similar to that used by PDH, and we refer 
the reader to that paper for a discussion of many of the 
assumptions involved. As a starting point, we also assume 
the same qualitative shock structure in the WBB as shown 
in Fig. 1 of PDH: that is a central wind source surrounded 
by a region of unshocked wind, separated from an outer 
region of shocked wind by an inner shock. The swept-up 
ambient medium is assumed to radiate efficiently and col- 
lapse to a negligibly thin shell coincident with the contact 
discontinuity sepa ratin g the stellar and ambient gas (cf. 
Dyson & de Vries |l972|) . 

A central and fundamental difference between the con- 
ductive case considered in PDH and the ablative case ex- 
amined here, however, is that here mass loading may also 
occur in the unshocked stellar wind. The effect is to heat 
this part of the flow, which leads to reduced Mach num- 
bers, and weaker jump conditions across the inner shocks. 
We shall see that for high mass loading rates, it is possible 
for the unshocked mass-loaded wind to connect directly 
with the contact discontinuity, without the presence of an 
inner shock. 

In this work we assume that the flow can be treated 
as a single fluid. This requires that the ablated clump ma- 
terial merges with the global flow in the sense that its 
temperature, velocity, and density approach those of the 
surrounding tenuous material. Ablation by itself might re- 
quire that the clump material mix microscopically with 
the wind to reach the density and temperature of the 
wind. However, though the mass-loss rate may be con- 
trolled by ablation, thermal conductivity will almost cer- 
tainly be responsible for material, once it is stripped from 
a clump, reaching the physical state of the tenuous fast 
flowing medium. Thermal conduction accomplishes this 
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phase transition without microscopic mixing, and acceler- 
ation to the global flow speed is effected by the response 
of stripped material to pressure gradients and viscous cou- 
pling, which may arise from a host of mechanisms includ- 
ing turbulence. Thus we envisage a two-stage process in 
which ablation controls the rate at which mass is stripped 
from the clump but conductivity becomes important for 
the merging into the global flow. 

This receives support from the comparison of the con- 
ductively driven evaporation time of a spherical clump 
with nT = 10 7 cm" 3 and T = 10 4 K, embedded in a 
medium with the same pressure (which is typical of plan- 
etary nebulae, Wolf-Rayet wind-blown bubbles, and star- 
burst superwinds) and T — 10 7 K (which is also typical of 
hot material in such regions) to the sound crossing time in 
the clump. The sound crossing time is somewhat (but not 
hugely) smaller than the conductively driven evaporation 
t ime f or a large range of clump sizes, if the Cowie & McKee 
( 1977 ) estimate of the mass- loss rate driven by saturated 
conduction is used. Consequently, ablation initiates mass- 
loss because it causes an increase in the surface area of a 
clump triggering more conductive heat transfer. Use of the 
analysis of Cowie & McKee (1977) and McKee & Cowie 
(1977) shows that for the assumptions given above, radia- 
tive losses do not hinder conductively driven evaporation 
unless the clump radius is greater than rs 15/n/ pc where 
rif is the number density (cm -3 ) of the global tenuous 
flow. Of course, a clump that does not cool radiatively af- 
ter it is compressed by a shock will have a shorter sound 
crossing time and be ablated more rapidly. However, by 
assuming a clump temperature of 10 4 K above we have 
established that the two-stage process is likely to be rel- 
evant in many environments if a clump radiatively cools 
after being shocked. 

Observations of WBBs and PNe indicate that mass 
loading may not necessarily begin at zero distance from 
the wind source. For instance, in the young nebula PN 
Abell 30 the clumps appear alm ost all of the way down to 
the central star (Borkowski et al. 1995 ), whilst in the more 
evolved Helix Nebula the clu mps a re all located closer to 
its edge (e.g. Meaburn et al. 1996 ). We therefore include 
the radius at which mass-loading 'switches on' as a free 
parameter in our models, and assume that the wind un- 
dergoes free expansion interior to this radius. Hence each 
solution will consist of a region of supersonic wind with no 
mass loading and an adjacent region with mass loading. 

One can imagine two possible causes for this minimum 
'mass-loading' radius. In one scenario the clumps could 
have been ejected at low velocity from the central star at 
an earlier evolutionary stage. The ejection of clumps then 
abruptly stopped, so that at the time of observation they 
had travelled a finite distance from the central star. By 
this process a central region clear of clumps surrounded 
by a clumpy region can be generated. A second possibil- 
ity is that clumps interior to the mass-loading radius have 
been completely destroyed by the wind. It seems reason- 
able to suppose that clumps located closest to the central 
star will be destroyed first, since they will have been sub- 



jected to the wind from the central star for the longest 
time. Then as the bubble or nebula evolves, clumps at 
ever increasing distance from the central star will be de- 
stroyed. Timescales for the destruction of clumps by ab- 



lation can be estimated from Hartquist et al. (1986) and 
Klein et al. (1994). Estimated destruction timescales vary 
from significantly less than to greater than the age of the 
bubble/PNe, in accord with the different spatial distribu- 
tion of clumps in objects of differing age. 

Regardless of which of the above scenarios is respon- 
sible for the existence of such a minimum mass-loading 
radius, this radius will physically increase with time. Our 
similarity solution requires that it increases in the same 
way as that of the contact discontinuity (i.e. r ex f 2 /( 5+A ), 
where A is the radial dependence of the mass- loading). For 
most of the solutions presented in this paper, the mini- 
mum mass-loading radius scales with or close to t. Since 
on physical grounds we might expect it to scale as t, our 
solutions closely match this requirement. 

In our solutions an inner shock may or may not be 
present - in the latter case the mass-loaded wind directly 
connects to the contact discontinuity, and the mass load- 
ing may be strong enough for the wind to become subsonic 
with respect to the clumps before the contact discontinuity 
is reached. If an inner shock is present, the postshock flow 
is by definition subsonic with respect to the shock, but 
may still be supersonic with respect to the clumps. In this 
case a number of different profiles for the Mach number 
are possible before the contact discontinuity is reached. 

At the center of the bubble prior to the onset of mass 
loading, we solve only the continuity and momentum equa- 
tions, with the implicit assumption that the thermal en- 
ergy of the flow is negligible, whilst in the mass loading 
regions we additionally solve the energy equation, and in- 
clude a source term for mass injection in the continuity 
equation. For a 7 = 5/3 gas, the equations for the mass- 
loaded flow are: 



dp 

at 



+ v • (H = p 



d{pu) 



dP 



dt 

d{\pu 2 + e) 
dt 



V • (pu') + ^-=0 



Or 



+ V-(^(^ 2 + |-)) = 
2 dp 



(1) 
(2) 

(3) 



In these equations the symbols have their usual meanings. 
In the next sections we discuss appropriate similarity vari- 
ables for these equations, our treatment of the boundary 
conditions, and the scaling relationships to normalize the 
resulting solutions. The reader is again referred to PDH 
for a more in-depth discussion of the details. 

2.1. The Similarity Variables 

Let the interclump ambient medium have a density of 
the form p = p^r 13 , and let us consider the case in 
which the mass-ablation rate is also radially dependent: 
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p = Qr x M 4 / 3 for subsonic ablation (M < 1), and p = Qr x previously discussed, this can sometimes occur before the 

for supersonic ablation (M > 1) (cf. Hartquist et al. 1986Q . specified position of the inner shock (x — Xi s ) is reached, 

A similarity solution demands that A = (2/3 — 5)/3, and in which case the mass-loaded flow connects directly to the 

the 'physical' parameters r, p, u and e may be expressed CD. If an inner shock is present, its velocity with respect 



in terms of the dimensionless similarity variables x, f(x), 
g[x), and h(x) where 



r = s Q-l/(B+A)£l/(6+A) t 2/(B+A) 

p = Q 5 /(5+A)£;A/(5+A) t (5+3A)/(5+A) j,s 
U = Q-V(5+A)£;l/(5+A) t -(3+A)/(5+A) g(a , ) 
e = Q3/(5+A )i; (2 + A)/(5+A) t -(l-A)/(5+A) ft(:c) _ 



(4) 
(5) 
(6) 

(7) 



Upon substituting the above similarity variables, the hy- 
drodynamic equations for the region of freely expanding 
wind become: 



0-^)/' 

(.9 - 5Tl) f + (2/ 



f9' + HF 

2xf 
(5+X)g 



5+3A 
5+A 



/ = o 



2/£ , 2+2A f 
x ' 5+A J 



(8) 
(9) 



where a prime denotes derivation with respect to x. For 
the flow in the mass-loaded region we obtain: 



to the bubble center is g s — 2xi S /(5 + A), and the Mach 
number of the preshock flow with respect to the shock is 



Mi = (,9i -g a )\ 



h 



7(7-1)^1 



(13) 



where the subscript '1' indicates the preshock values. 
Because (unlike for the evaporative case) M\ is not neces- 
sarily large, we cannot assume that the strong shock jump 
conditions apply and therefore must use the full Rankine- 
Hugoniot relations. The postshock values are then: 



h = 



7' 



1 M 2 



7 - 1 M 2 + 2 



h 



92 = Xis 



(51 -xi S )-r 
h 



2 7 M 2 - (7 - I] 

7+1 



hi. 



(14) 
(15) 
(16) 



3A 



( 9 - ft) /< + /.,< + ^ + ^/ 



x 5 + A 

A M 4 / 3 = 



(■9-5+t) 



{9 - 5+t) 



3 ? 5 + A' 



h' 



hg' - 



9 A M 4/3 
1 



(10) 



(11) 



A. 



5 + A 
10% 
3 x 



9 —x x M^ 3 = (12) 



where the Mach number, M — 3\/9//(10/i). It is simple 
to rearrange these equations to find /', g' , and h' which 
may then be integrated to obtain solutions. 

2.1.1. Boundary Conditions 

The values of / and g at x = Ax, the inner boundary 
of our integration, are g = \/2/<j) and / = l/(27r Ax 2 g 3 ), 
where is a parameter which specifies the relative amount 
of mass loading in the solution, h is set to zero at this 
point. Equations N and [| are integrated to the radius at 
which mass loading begins. From this point onwards we 
then integrate Equations [HH13 instead. The Mach num- 
ber of the wind rapidly decreases from its initial value of 
infinity and both the density and velocity of the wind re- 
spond to the mass addition. Integration proceeds until the 
contact discontinuity (CD) is reached. This occurs when 
the flow velocity is equal to the coordinate velocity (i.e. 
v = dR/dt), which is where g{x = x c d) = 2x/(5 + A). As 



Even if the flow velocity has not converged to the coor- 
dinate velocity, an inner shock may not exist. This can 
occur if the mass loading interior to the proposed position 
of the inner shock is strong, resulting in Mi being less 
than unity. In such a scenario, the jump conditions are 
not applied and the integration proceeds as before. 

At the CD, we must also satisfy conservation of mo- 
mentum. We take 



Q-l/(5+A)£l/(5+A) 



P(} 



1/(5+0) £i/( 5+/3) 



(17) 



as a measure of the ablated mass to the mass in the swept- 
up shell. It follows that at the contact discontinuity, 



V(5+/3) 







,1 - 1)' 1 


T fl 


g 2 - 






(5+A)(3+/3) X yJ 



(18) 



The rate of massloss from the star is 



M c = lim Am x fgz, 

x—>0 



where 

S = Q 2/(5+A) E {3+X)/i5+X) i(S+2A)/(5+A) j 



(19) 



(20) 



which differs from the definition in the conductively driven 
case. The mass loading parameter <f>, which is a measure 
of the ratio of ablated mass to wind mass is given by <f> — 
M c /E. 
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Table 1. The influence of x m \ on solutions with A = —3. Two cases are considered: minor mass loading (0 = 100) 
whereby the value of x m i is relatively unimportant, and appreciable mass loading (0 = 10) where x m i exerts a major 
influence. From left to right, the columns indicate: i) the value of x rn i, ii) the ratio of the shock radii and CD, iii-vi) 
the fractions of energy that are thermal and kinetic in the bubble, kinetic in the swept-up shell, and radiated, vii) 
the value of 9, viii) the ratio of mass within the bubble to the total mass injected by the wind (values greater than 
unity indicate mass loading), ix) the ratio of mass in the swept-up shell to mass in the bubble, x) the preshock Mach 
number of the unshocked wind. 
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2.2. Scale Transformation and Normalization 

The dimensionless Eqs. ^0[]l^ are invariant under the fol- 
lowing transformation, which we shall call a normaliza- 
tion: 



(21) 



These can be combined to obtain the normalizations for 
the mass, and for the kinetic and internal energies of the 
bubble: 



x — ► 


OX 


/-» 


a x f 


9-> 


ag 
a 2+x h 



m = Ait I fx dx 
k = 4n -fg x dx 
4ir / hx dx 



,3+Ai 



,5+Ai 



o 



5+Ai 



The full equation for the mass in the bubble is: 



,2 J„ _ „3+\ IVI rl 



M b = 4tt a i+A Et fx z dx = a A+A —-m. 



x=0 



The total mass lost by the star during its lifetime is 

5 + A 



AL 



w 



M{t)t 



11 + 3A 



-Met. 



(22) 
(23) 
(24) 

(25) 
(26) 



The degree of mass loading in the bubble, ${,, is defined as 
the ratio Mb/AL W . The mass-swept up into the shell, M s h, 
the ratio of swept-up mass to bubble mass, M s h/ALt, the 
kinetic energy of the shell, KE s h, and the pdV work on the 
shocked ambient gas, are all identical to the corresponding 
forms in PDH. 



2.3. Solution Procedure 

If an inner shock exists, three main parameters determine 
the form of a given similarity solution (A, and Xi S /x c d). 
However, the ratio of the radius of the onset of mass load- 
ing to the radius of the inner shock is an additional param- 
eter unique to the ablative solutions, and its value (x m i) 
may also influence the form of the solution. We find that 
if the mass loading of the bubble is minor, x m i has rel- 
atively little influence, and vice-versa (since p oc r A and 
~~ 3 < A < -2 for most of the solutions presented, most 
of the mass-loading occurs close to the minimum mass- 
loading radius). Alternatively, if an inner shock does not 
exist, Xi s /x c d and x m i have no meaning as we have defined 
them. In this case we report the ratio of the onset radius 
of mass loading to the radius of the contact discontinuity 
as the parameter x m \. 

The ratio of mass pick-up from the clumps to the mass 
swept-up from the interclump medium (which is related to 
the value of 9) is obtained only after a particular solution 
has been found. For bubbles whose evolution is signifi- 
cantly altered by mass loading, we require that simulta- 
neously and 9 are both small. 

The similarity equations were integrated with a fifth- 
order accurate adaptive step-size Bulirsch-Stoer method 
using polynomial extrapolation to infinitesimal step size. 
Once the CD was reached, the similarity variables were 
rescaled using the relationships defined in Eq. |2l] so that 
x c d = 1. The mass, and kinetic and thermal energies of 
the bubble were calculated, as were the kinetic energy of 
the shell and the energy radiated from it. The correct nor- 
malization to satisfy global energy conservation was then 
obtained. Finally, for given values of E, Q and t, the sim- 
ilarity variables x, /, g and h may be scaled into the phys- 
ical variables r, p, u and e. 
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Fig. 1. Results for A = —3, Xis/^cd = 0.687, and (/> = 50, 000. In the limit of large 4> (i- e - negligible mass loading), the 
earlier results of Dyson (1973) are recovered. From the top the panels indicate density, velocity, internal energy, and 
Mach number (with respect to the stationary clumps) as a function of radius. Except for the Mach number, the panels 
are normalized to values of 1.0 at the contact discontinuity. 



3. Results 



We first checked our work against the solutions obtained 
by Dyson (1973) for WBBs with no mass loading. For 
this comparison it was required that A = — 3 and that 
4> was large. For A = — 3 the radius of minimum mass 
l oadin g increases as t. Excellent agreement with Dyson 
( 1973 ) was found over a large range of Xi S /x c d- For large <f> 
(i.e. negligible mass loading) the value of x m \ has no effect 
on the resulting solution. Fig. [y shows a sample solution. 



In Fig. H we show results for A = —3, <f> = 100, with 
different values of x m i . The mass loading in these results 
is small, but not negligible, so that the precise value of 
x m i has minor consequences for the solution obtained. In 
Table H we tabulate important parameters from these so- 
lutions. In particular, we find that values for the ratio 
Xi s /x c d and the energy fractions are very similar. The frac- 
tional mass loading of the bubble, <!>;,, is somewhat more 
sensitive to the value of x m i, and increases when the 'turn- 
on' radius for mass loading is decreased, as one would ex- 
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Fig. 2. Results for A = —3, <fi = 100, and various values of x m i- For these parameters the mass loading is small but 
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Fig. 1. In Fig. 3 we show solutions where the value of x m i has a much greater influence. 



pect. The preshock Mach number also reflects the degree 
of preshock mass loading, again as expected. Interestingly, 
the profiles of postshock Mach number are, however, very 
insensitive to the preshock Mach number. With the as- 
sumption of constant <j>, we note that since Xi S /x c( i varies 
for the above results it was important to see whether this 
was a possible factor for some of the differences in these 
solutions. Therefore we also investigated the effect of dif- 
ferent values of x m i whilst keeping Xi S /x c d constant and 
varying <f> as necessary. We find that the resulting solu- 



tions are again fairly similar, and differ at about the same 
level as the results in the top half of Table H. Therefore 
it seems that varying x m i has a direct effect on the re- 
sults, and not just an indirect influence through changes 
induced in Xi S /x c d and/or <f>. In the following, therefore, 
we will vary x m i whilst keeping <j> fixed. 

At the other extreme, Fig. y shows results where the 
value of x m i can have large consequences for the solutions 
obtained. This occurs for values of <j> for which the result- 
ing mass loading of the bubble is important. Parameters 
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from these solutions are tabulated in the lower half of 
Table |lj. Although the ratio Xi S /x c d remains relatively in- 
sensitive to the value of x m i, the fractional mass loading 
$b, the ratio of the shell mass to the bubble mass M s h/Mb, 
and the energy partition are all significantly affected, as 
Table [y shows. This is expected given that most of the 
mass loading occurs near the minimum mass loading ra- 
dius for A = —3. 

An interesting consequence of mass loading the wind 
prior to the innershock is that if the mass loading is large, 



the Mach number of the flow can be reduced so much 
that the flow directly connects to the contact discontinuity 
without the presence of an inner shock. In Fig. |l| we show 
one solution where this happens. In this example, the flow 
is continuously supersonic with respect to the clumps, and 
we obtain $b = 1.25 and M s h/Mf, = 0.07. This perhaps 
indicates that this morphology is most favoured at early 
stages of the evolution of the bubble before a significant 
amount of ambient medium is swept-up. In another ex- 
ample of 'direct connection', we find a solution where the 
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Fig. 4. Results for A = —3, <f> = 50. Here mass loading is so severe that the flow directly connects to the contact 
discontinuity. For the solution shown, the onset of mass loading occurs at x = 0.359iEcd- Once the flow starts mass 
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Mach number of the flow with respect to the clumps drops 
below unity (see panel f of Fig. |g). In this case, 4>b = 3.4 
and Mgh/Mb — 480, suggesting a much older bubble. We 
also find that low values of A (~ —3) are apparently needed 
to obtain directly-connected solutions. 

In Fig. M we compare solution profiles for the case 
where A = — 2 and for negligible and high mass loading 
rates. A A = — 2 clump distribution is the most phys- 
ically plausible for many astrophysical objects. For this 
distribution the minimum mass loading radius increases 
as t 2 ' 3 , which is not far removed from a linear increase 



with t. The solutions have been scaled from the similarity 
results by assuming current values for the wind param- 
eters appropriate for a WR star: M c = lO~ 5 M0yr -1 
and v w — 1500 km s -1 . The mass- loss rate increases with 
time as t 2 ^ 3 whilst the wind velocity decreases as i -1 / 3 , 
so the wind has evolved from a case more suitable to an 
O-star. The assumed bubble age is t = I0 4 yr. The in- 
terclump medium varies as p ex r -1 ' 2 (i.e. (3 = —1/2) 
which is a flatter distribution than for a constant velocity 
wind (r~ 2 ). However, since the wind mass- loss rate is in- 
creasing and the velocity decreasing, an r -1 ' 2 interclump 
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density profile is physically realistic for such a scenario. 
In comparison with the negligibly mass-loaded solution, 
the increase in density inside the bubble for the highly 
mass-loaded solution can be clearly seen, along with a de- 
crease in preshock velocity, preshock Mach number, and 
postshock temperature. Interestingly, the Mach number of 
the postshock flow is higher for the heavily mass-loaded 
solution. Note also that the heavily mass-loaded solution 
extends to a larger bubble radius. This is a consequence of 



the different preshock interclump densities: for the heav- 
ily mass-loaded solution n e = 0.066, whilst n e = 0.195 for 
the negligibly mass-loaded solution. 

In PDH it was shown that there is a negative feed-back 
mechanism caused by the evaporation of mass from em- 
bedded clumps, which set a maximum limit to the amount 
that a bubble could be massloadcd. Additionally, for the 
bubble mass to be greater than the mass of the swept-up 
shell, a large radial dependence of the mass loading was 



Pittard et al.: Self-similar evolution of mass- loaded WBBs 



11 




Fig. 6. The Mach number as a function of x for a number of different solutions. With respect to the clumps, the shocked 
region can be either entirely subsonic, entirely supersonic, or have one (or maybe more) sonic points. The panels show: 
a) shocked region entirely subsonic with a monotonic decrease in Mach number (A = —2, <f) — 0.5, Xi S /x c d — 0.585, 
x m i = 0.1); b) shocked region entirely supersonic with a monotonic increase in Mach number (A = —3, <f> = 100, 
Xis/x c d — 0.902, x m i — 0.9); c) postshock Mach number initially subsonic and decreasing, before levelling off and 
subsequently increasing through a sonic point (A = 4, <j> = 10, Xi S /x c d — 0.297, x m i — 0.1); d) shocked region entirely 
subsonic with a decreasing Mach number followed by a modest increase close to the contact discontinuity (A = —3, 
<f> = 10 5 , x is /x cl i = 0.181, x m i = 0.1); e) postshock flow continuously supersonic though with an initial decrease 
in Mach number (A = —3, <f> — 10, Xi S /x c d — 0.816, x m i — 0.02); f) continuously supersonic flow (solid) directly 
connected to the contact discontinuity (i.e. no inner shock) (A = —3, = 50, x m i = 0.359) and a flow becoming 
subsonic (dots) directly connected to the contact discontinuity (A = —3.2, = 12.5, x m i — 0.246). 



required (A > 4). In this work we instead find that small 
values of A are required to satisfy this condition. Our sim- 
ulations show that it is impossible to obtain a similarity 
solution with M s h < Mb for A > — 1 (e.g. for A = 1, 
M s h ^ 2.4Mb). However, for smaller values of A, solutions 
satisfying M s h < Mb can be found (e.g. for A = —2 (-3), 
we can obtain M b » 8 (29) M sh , where $ fc w 7.6 (15)). For 
a given value of M s h/Mb there appears to be a maximum 
value for ${,. It occurs when mass loading starts very close 
to the bubble centre (i.e. x m i — > 0) and when the flow re- 
mains supersonic relative to the clumps over the entire 
bubble radius. 

We have also investigated the Mach number profile of 
the solutions that we obtain, of which some are already 
shown in Figs. @ and |l As for the evaporative mass loading 
in PDH, we note that a number of different forms are 
possible, and summarize these in Fig. pi 

Finally, radial profiles of the X-ray emissivity per unit 
volume and temperature have been calculated (Fig. 0). 
We assume that the emissivity A oc 7i 2 T~ 1 / 2 , which is a 



good approximation over the temperature range 5 10 5 K 



< 



T < 5 10 7 K (cf. Kahn 1976). The upper panels show re- 



sults where an inner shock is present, and we will discuss 
these first. The solution with the solid line displays rela- 
tively constant values of A and T, and is characterized by 
a high degree of mass loading ($h = 9.7). Note that the 
highest temperature in the bubble does not occur post- 
shock, but rather shortly after mass loading begins. The 
very high emissivity at the onset of mass loading is a con- 
sequence of the low temperature at this point and the 
assumed T~~ 1 ' 2 dependence of our emissivity, and is not 
physical. The solution with the dotted line has an inte- 
rior dominated by free-expansion of the wind source. Mass 
loading only occurs over the final 20 per cent of the radius, 
and is relatively weak, and the maximum temperature in 
the bubble occurs immediately postshock. The solution 
with the dashed line has some characteristics of each of 
the other solutions. It is mass- loaded from small radius 
which again leads to the maximum bubble temperature 
occurring near the center. However, the mass loading is 
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Fig. 7. Radial profiles of emissivity per unit volume and temperature, normalized to values of 1.0 at the limb. The 
upper panels show results where an inner shock is present, whilst the bottom panels show results where the mass 
loading region directly connects to the contact discontinuity. For the upper panels the results shown are for: A = —2, 
</> = 0.5, Xi S /x c d = 0.59, x m i = 0.1 (solid); A = —3, <p — 100, Xi S /x c d = 0.90, x m i — 0.9 (dots); A = —3, 4> = 10, 
Xis/x c d = 0.82, x m i = 0.02 (dashes). For the lower panels the results are for: A = —3, <fi = 50, x m i = 0.36 (solid); 
A = -3.2, <j) = 12.5, x m i = 0.25 (dots). 



less severe (${, = 6.5) and its radial dependence steeper 
(A = —3) than for the solution with the solid line, so that 
there is a large fall in density between the switch-on ra- 
dius for mass loading and Xi S . This leads to a steep rise in 
emissivity from the inner shock to the bubble center in a 
similar fashion to the solution with the dotted line. 

The bottom panels of Fig. M show results where there 
is no inner shock and where the mass loading region di- 
rectly connects to the contact discontinuity. The Mach 
number profiles for these solutions are displayed in the 
bottom-right panel of Fig. ||. As can be seen, the temper- 
ature profiles are approximately flat over their respective 
mass loading regions, and are very similar in shape to each 
other. This contrasts with the behaviour of their emissiv- 
ity profiles, where it is clear that the solution represented 



by the solid (dotted) line has A generally decreasing (in- 
creasing) with radius. 

A general result from Fig. |7|, which was also found in 
the evaporatively mass-loaded bubbles of PDH, is that 
higher mass loading tends to decrease the central emis- 
sivity and increase the central temperature relative to the 
limb. 



3.1. Comparison with observations and other 
theoretical work 

Although our work cannot be directly com pared to th e 
hydrodynamic simulations of Arthur et al. ( 1993 , 1996 ), 
a number of similarities exist. These include a sharp in- 
crease in temperature at the onset of mass loading, and 
often dp/dx > over the mass-loaded region. A notice- 
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able achievement of the work of Arthur et al. (1996) was 
the reproduction of the observed correlation between the 
velocity and ionization potential of ultraviolet absorption 
features towards the central star (Smith et al. 1984). Our 
solutions also often reproduce this correlation (e.g. see 
Fig. a which shows the velocity and temperature decreas- 
ing together for the heavily mass-loaded solution), as do 
those of Hanami & Sakashita (1987). Interestingly, solu- 
tions with mass-loading from conductive evaporation of 
the swept-up shell do not reproduce this correlation (see 
Weaver et al.|1977|). 



temperature. The maximum temperature in the bub- 
ble is often not the post-inner-shock temperature, but 
occurs near the onset of mass loading. In some cases 
this can be many times greater than the post-inner- 
shock temperature. 
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4. Summary 

We have investigated the evolution of a mass-loaded wind 
blown bubble with a constant rate of injection of mechan- 
ical energy from a central wind source. The mass loading 
occurs due to the hydrodynamic ablation of distributed 
clumps, and is 'switched-on' at a specified radius, interior 
to which the wind expands freely. The requirement that 
the solution be self-similar imposes a link between the ra- 
dial variation of the interclump density (p oc r ) and the 
rate of mass loading from the clumps (p oc r A ) which forces 
A = (2/3-5)/3. 

We first produced solutions with negligible mass load- 
ing and /3 = — 2 (which correspond to constant M W i n d 
and v W j nd), w hich we compared with results obtained by 
Dyson ( 1973| ). Excellent agreement for the structure of 
the bubble was found. We also confirmed that for negli- 
gible mass loading, the value of x m i had no effect on the 
resulting solutions, as desired. 

We then investigated the changes in the structure of 
the bubble for different values of A, (j), X[ s /x c d, and x m i. 
The central conclusions are: 

— Substantial mass loading of the wind-blown bubble can 
occur over a wide range of A. However, to additionally 
satisfy the requirement that the bubble mass is larger 
than the mass of the swept-up shell, low values of A 
are needed (e.g. A < —2). 

— The profiles of the flow variables are significantly al- 
tered under conditions of large mass loading. With 
respect to solutions with negligible mass loading, the 
density and temperature profiles increase, and the ve- 
locity and Mach number profiles drop. Changes to the 
profiles can be very rapid. Mass loading of the wind 
also reduces the Mach number prior to the inner shock. 

— The mass-loaded wind may also connect directly to 
the contact discontinuity without the need for a global 
inner shock. 

— The Mach number of the flow relative to the clumps 
which are injecting the mass can take several different 
forms. The flow can be either entirely supersonic, or 
have one (or maybe more) sonic points. 

— Whilst mass loading reduces the bubble temperature, 
it also tends to reduce the emissivity in the interior of 
the bubble relative to its limb, whilst simultaneously 
increasing the central temperature relative to the limb 
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